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ABSTRACT 


Final Technical Report on NASA Grant NSG 1375 (January 1977 - March 1978), 
Processing Air Quality Data : The Virginia State Air Pollution Control 

Board (VSAPCB) has conducted air quality sampling since January 1968. Data 
for air quality parameters, recorded on written forms, punched cards or 
magnetic tape, are available from the VSAPCB Region VI ~ Hampton Roads for 
1972 through 1975. The principal investigator for this research project 
participated at the Langley Research Center (LRC' as a National Aeronau- 
tics and Space Administration - American Society for Engineering Education 
(NASA-ASEE) 1976 Summer Research Fellow in the initial phase of computer 
software development by the LRC Analysis and Computation Division (ACD) 
for the processing of air quality data. During this initial phase, com- 
puter software was developed to (1) calculate several daily statistical 
measures of location, (2) plot time histories of data or the calculated 
daily statistics, (3) calculate simple correlation coefficients, and (4) 
plot scatter diagrams. This research project continued at Mississippi 
State University (MSU) the development of computer software for process- 
ing air quality data to include time series analysis and goodness of fit 
tests. The usual multiplicative time series model with trend, seasonal, 
cyclical and irregular components, Y = TxSxCxI, was extended to include 
daily and weekly periodic components, Y = TxDxWxSxCxl. Computer soft- 
ware was developed to (1) calculate a larger number of daily statistical 
measures of location, (2) calculate a large number of daily, monthly and 
yearly measures of location, dispersion, skewness and kurtosis, (3) decom- 
pose the extended time series model and (4) perform some goodness of fit 
tests. The computer program is described, documented and illustrated by 
examples. Recommendations are made for continuation of the developement 
of research on processing air quality data. 
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I. BACKGROUND 


Air quality data has been collected by the Virginia State Air Pollu- 
tion Control Board (VSAPCB) since January 1968. Among the seven VSAPCB 
regions, Region VI - Hampton Roads encompasses the Norfolk-Peninsula area 
from the North Carolina border to the Williamsburg-Yorktown vicinity. 

The National Aeronautics and Space Administration (NASA) has entered into 
a cooperative agreement with the VSAPCB, by and through tie Region VI 
Director, whereby the Langley Research Center (LRC) assists the Region VI 
office by making available the technological capabilities of the LRC that 
may be of value to the Region VI office in accomplishing the VSAPCB' s 
monitoring and emissions control programs. 

In addition to the considerable expertise of the LRC for assisting 
the Region VI office of the VSAPCB in the realm of sampling instrumenta- 
tion, the vast computer complex, computer programming, and mathematical 
analysis capabilities of the LRC Analysis and Computation Division (ACD) 
can be of invaluable assistance in the data processing of very large 
amounts of air quality data collected and stored over the past several 
years. 

The principal investigator (Dr. Charles M. Butler) for this research 
project had the privilege of participating as a NASA-ASEE (American Society 
for Engineering Education) Research Fellow in the initial phase of the 
LRC effort to assist in the analysis and interpretation of Region Vi's 
accumulated air quality data. During this initial phase (Summer 1976), 
computer software was developed to access the stored data, to calculate 
some summary statistics and correlation coefficients, and to plot some 
time histories and scatter diagrams. 


The proposed purpose of this research project was to continue the 
development of computer software for analysis and interpretation of air 
quality data. 

II. RESEARCH PROPOSED 

A. INTRODUCTION 

Any complete description of data processing must include the five 
steps : 

1. generation, 

2. recording, 

3. validation, 

4. storage and retrieval, and 

5. analysis and interpretation^. 

A great deal of effort has been and continues to be expended on the first, 
second and fourth of these steps. Step three probably needs careful 
study but was not the subject of this research project. Step five, the 
analysis and interpretation of air quality data, was the topic to be 
investigated herewith. 

During the summer of 1976, the principal investigator of this pro- 
posed research project participated in the initial phase of the LRC-VSAPCB 
effort to develop software for meaningful analysis and interpretation of 
air quality data collected in Region VI over the past several years. This 
data consists of continuous hourly averages of sulphur dioxide CSO 2 ) , 
carbon monoxide (CO), ozone CO^) , and wind velocity; continuous two-hourly 
averages of soiling index (coefficient of haze - COH) ; and daily averages 

^Superscripts refer to numbered references in Section V of this report 



of total suspended particulate (TSP) , collected either every third day 

or every sixth day. Such data may be recorded on written forms, punched 

cards, or magnetic tape; and it may be in standard or non-standard formats. 

1 

Since all such data may be transformed to the standard SAROAD formats 
and may be recorded on magnetic tape, the software was to be developed 
for accessing data on magnetic tape in the standard SAROAD formats* 

As of August 13, 1976, under the supervision of Dr, Stephen K, Park 
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of leg’s ACD, a software package had been developed at LRC to plot (1) 
continuous hourly or two-hourly time histories for each month; (2) daily 
arithmetic mean, minimum, maximum, midrange and median for each month; 
and (3) correlation scatter diagrams for pairs of parameters at each site 
and at pairs of sites for each parameter for each month. 

This research project was proposed to continue development of this 
software package to (1) decompose the data for each parameter into long 
term trend, seasonal variation, cyclical fluctuation and irregular 
components and (2) perform goodness of fit tests for various standard 
distribution models to the air quality data. 

This research project was proposed to consist of two parts: time 
series analysis and goodness of fit tests, 

B. THE TIME SERIES ANALYSIS 

^ time series of the type considered here is a sequence of values 
of some variable (or composite of variables) arranged in the chronological 
order of their occurrence. Each value (datum) is assumed to consist of 
four components : 

Trend, denoted by T, is the long-term (many years) smooth variation 
over time. Most often T is considered linear, but is of other forms if 
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the situation or date so indicates; polynomials, exponentials, and power 
functions are not uncorauon. The method of least squares is usually the 
method used to fit the trend curve to the data, although experienced 
practitioners may frequently use visual or other short-cut methods for 
data with which they are thoroughly familiar. 

Seasonal variations , denoted by S, is a conglomerate of all periodic 
variations in the data whose periods are one year or less. Usually S is 
considered in the form of an index number which expresses seasonal varia- 
tions as a percentage of the corresponding trend (or combination of trend 
and other components) values. There are two basic methods of extracting 
the seasonal variations component from time series data of the type con- 
sidered here: the ratio-to-trend and the ratio-to-moving-average methods. 
I'/here there is not yet enough experience with a situation or type of data 
to know which method produces the "best" results, both methods should be 
applied and their results compared to see which produces the better 
representation of the data or situation being studied. 

Cyclical fluctuations , denoted by C, is a conglomerate of all peri- 
odic fluctuations in the data whose periods are greater than one year — in 
practice, usually up to about three to seven years, depending on the nature 
of the data. C is usually expressed as an index number giving the cyclical 
component as a percentage of trend or combination of trend and seasonal 
variations. The ratio-to-moving-average method is used to extract the 
cyclical component in situations where it is believed to be of sufficient 
importance, but, quite often, it is not considered as a separate entity, 
being left in combination with S when using the ratio-to-trend method. 

Irregular or random movements , denoted by I, consists of movements 
in data values that are not attributable to trend, seasonal variations 


or cyclical fluctuations* There is usually no effort made to analyze the 
irregular component, but only to recognize it where it occurs in the data 
and remove or minimize its effect on the data values. Often the I com- 
ponent is attributable to tornado, hurricane, lightning, flood, fire or 
other rare and unpredictable phenomena. 

The time aeries model consists of the combination of the four com- 
ponents in some multiplicative or additive fashion, although various 
combinations and transformations may be appropriate in some unusual cir- 
cumstances. The model most often used is the straightforward multiplica- 
tive model Y = TxSxCxI, where each datum is considered to be the simple 
product of the four components. 

If the ratio- to- trend method is used, a trend line is fitted to the 
data, usually by the method of least squares. The trend value correspond- 
ing to each datum y-value may then be calculated by substituting that 

A 

datum x-value into the expression for the trend line, T = a + bX 
(assuming a linear trend — if other than a linear trend line is used, 
the same procedure is followed with that fitted trend curve). Each trend 
value is then divided into its corresponding datum y-value (In terms of 
the model this is = SxCxI.) yielding an index number (It is customary 
to multiply by 100 to discuss the index number in terms of percentages.) 
which expresses the SxCxI component of each datum y-value as a pro- 
portion of the corresponding trend value. In the ratio- to-trend method 
this quotient is referred to as the seasonal index ; this is rather "loose” 
terminology, since it ac' ally is an index of combined seasonal, cyclical 
and irregular movement. If there is sufficient data (more than one year) 
for repeated calculations of the same index (e.g., for March 1973, March 
1974, and March 1975), the repeated indexes are usually averaged to 
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produce one index for that point (month, day, hour, week, etc.) in the 
time series. The ratio- to- trend method then uses the trend line (curve) 
to forecast or predict other values of the variable being studied, usually 
in the near future (Of course, the usual precautions about extrapolation 
are in order.); these forecast or predicted trend values are then corrected 
for seasonal variation (actually combined seasonal, cyclical, and irregu- 
lar variations when the ratio-to-trend method is used) to produce the 
value expected if and when the actual data is collected. (The usefulness 
of such a forecasting ability for air pollution data would be to cut back 
on emissions to avoid the forecast value when it is beyond allowable 
limits. ) 

I If the ratio-to-raoving-average method is used, the trend line may 

i be fitted (usually by the least squares method) either before or after 

the seasonal component is extracted — or even after extracting both 

■ seasonal and cyclical components; it is not clear what effect the order 

> of the operations has on the final decomposition of a time series and a 

' selection is usually made in a previously unexplored situation after 

comparing the results of the various orders in which the operations may 

be performed. The seasonal component is extracted by dividing each datum 

y-value by the annual average centered on that datum. (In tems of the 
Y 

model this Is - — - — =- = S.). If it is desired to extract the cyclical 
1 X U X t 

component, each datum y-value is divided by a three to seven (or more, 

depending on the nature of the data) years average centered on that datura 

Y 

(In terms of the model this is ■ = S x C.) followed by division by the 

J. X 1 

S value found In the preceding annual ratio-to-moving-average calculation 
S X C 

( — ~ — =5 C). Then a trend value is forecast or predicted, corrected for 
seasonal or seasonal-and-cyclical variation (The cyclical correction is 
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ignored if the preceding analysis and calculations have indicated it is 
of little importance.)* and used as appropriate in the analysis. (In 
analysis of air pollution data, this may prove to be a valid method of 
predicting parameter values that are expected to exceed allowable limits, 
allowing a reduction in emissions to prevent the undesirable excess.) 

C. GOODNESS OF FIT 

The two commonly used tests for deciding how well a theoretical 

distribution fits a set of data reduced to a frequency distribution are 
2 

the X (chi'-square) and Kolmogorov- Smirnov goodness of fit tests. Both 

are appropriate for the problem of fitting a theoretical distribution to 

2 

air pollution data of the type considered here, but only the x will be 

described. (The Kolmogorov-Smirnov test could be applied later to see 

2 

how well it agrees with the x test.) 

The data is first reduced to a frequency distribution using relative 

frequencies and a number of classes appropriate for the data being studied. 

(There is no concrete rule but a good rule of thumb is to use between 

five and fifteen classes — exceptions to the rule are not uncommon and 

only a trial-and-error approach can decide the number of classes to be 

used in a previously unexplored situation.) 

A theoretical distribution is selected as a possible source of the 

data being studied, e.g., for total suspended particulate, the literature 

4 

indicates a lognormal distribution has been used by most recent prac- 
titioners. (I believe a study should be made, mostly empirical, testing 
the goodness of fit of a variety of distributions — including some modifi- 
cation of the Poisson — to decide what theoretical distribution to use for 
each of the air pollution parameters studied.) The probability for each 
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of the frequency distribution class intervals is obtained for the selected 
distribution (e.g, , lognormal for total suspended particulate) using pa- 
rameter estimates calculated from the sample data (frequency distribution). 

For example, if Y = e is lognormal, X is normal with mean p and variance 
2 2 

a , and p and a are considered to be the parameters of Y also. Then the 

,2 2 2 

mean of Y is the variance of Y is ^ “ 1) • 

_ n _ c 

The parameters would be estimated by p = X = .f-ln Y Jn = X' = .S-f^-'ln 

1~1 i i~X 


n 


,-tv 2 


Y'./n and by ^ ~ .Z,(ln Y. - X)^/{n - 1) ^ .1., (In Y'. - X’)^/(n - 1), 

Xi x[ 

where Y, = e for i = 1,2,..., n are the n data values: or Y! = e for 

i = 1,2,.,., c are the c class marks and f^ for i = 1,2,.,,, c are the c 

class frequencies; and In is an abbreviation for natural logarithm. (Since 

two values have been calculated from the data to obtain the theoretical 

2 

frequencies, the number of degrees of freedom for the x test, will be n - 2.) 

A standard normal distribution table is then used to obtain the theoretical 

frequencies (probabilities) for the classes in the frequency distribution, 

e.g., if f > j = ’,2,..., c is the observed frequency of the jth class 

which has lower real limit X._, and upper real limit X., the corresponding 

expected theoretical frequency for that class xi^ill be given by f^. = 

P(Y._ < Y ^ YJ = P(Y. < < Y ) = P(ln Y < X ^ In Y.) = P{ (In Y. - 

J J J J 

p)/a < (X - p)/o .< (In Y^ - p)/a} = ${ (In - X)/s) - ${ln Y,_j_ - X)/s>, 

using p = X and = S^, where $ is the cumulative normal distribution 
function. (Hence, any standard normal distribution table can be used to 
calculate lognormal probabilities. Of course, we are assuming very large 
sample sizes since we a' treating (X - X)/s as though it were normally 
distributed. ) 

2 

A value of x (observed) is then calculated from 

X^ = .2 (f . - f .)^/f .. 

oj ej ej 

2 

This value is then compared with tabulated values of x,^ with v = n - 2 
degrees of freedom to determine at what level of significance the hypothesis 
that Y is a lognormal distribution may be rejected. 
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III. RESEARCH ACCOMPLISHKP 
A. THE EXTENDED TI^^E SERIES ANALYSIS MODEL 

After seeing a few plots of air pollution data as time series of 
hourly or txTO-hourly averages, it is possible to surmise periodic com- 
ponents which would not be discernible when using the usual decomposition 
of the multiplicative time series model Y = TxSxCxI. Hox^rever, it is 
not difficult to extend this model to incorporate additional multiplica- 
tive components x^ith periods of any desirable duration. 

The extended time series model consists of the usual four components, 
T, S, C, and I, as described in I.B. above, together x^lth the follox^ing 
txjo components: 

Daily variations , denoted by D, is a conglomerate of all periodic 
variations in the data xjhose periods are 24 hours or less. Like S and C, 

D x^;ill be expressed in the form of an index number. 

Weekly variations , denoted by W, is a conglomerate of ail periodic 
variations In the data x^hose periods are seven days or less. Like S, C, 
and D. will be expressed in the form of an index number. 

The extended time series model , as considered here, x^rill be expressed 
as the multiplicative model given by V = TxDxWxSxCxI, and the 
ratio-to-moving-average method will be used to extract first D, and 
then W components. 

To extract D, centered moving averages using 24 hour periods xvlll be 
calculated for all possible 24 hour periods. This, In effect, eliminates 
or ''smoothes out" all components x^ith periods of 24 hours or less, leaving 
each centered 24 hour average consisting of the components TxWxSxCxI, 
Thus, if we divide each of the original Y data values by its correspond- 
ing centered moving average, the result xatIII be 
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Y _ T X D X W3C S X C X I _ 

TxWxSxCxI “ TxWxSxCxI “ 

which then is the desired index number for the daily variation. We 

would then average the D values for corresponding time periods each day> 

A 

e.g,, midnight to 0200 for two-hourly averages, and use this average-D 
to correct the forecast trend value for that time period. We do like- 
wise for each time period each day. 

To extract W, centered moving averages using seven day periods will 
be calculated, v7hich eliminates or "smoothes out" all components with 

A 

periods of seven days or less. Now, when we divide each Y data value by 
its corresponding centered seven day moving average, we have 

Y _ TxDxWxSxCxl ^ 

TxSxCxI “ TxSxCxI 

This result, in turn is divided by the average of the D values during 
that seven day period to get 

A A 

D X W 

— — = W, 

D 

avg. 

wh,ioh is the desired index number. Again we average the W values for 
each day of the week, and use that average-W to correct the trend fore- 
cast values which have already been corrected for the average-D. 

We may continue this process for S and C, in turn. 
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ORlGiNAt PAGE !S 

qP pqq^ quality 


The following data in the standard SAROAD format 


i 


is used for this 


example; 


1480710001F022730101001120181092999999999999999999999999999999999999999999999999 

1480710001F022730102001120181092999999999999999999999999011000800160006000800140 

1480710001F022730103001120181092018001800160012000400030006000800060007000400030 

1480710001F022730104001120181092003000300030003000600060006000500050005000609999 

1480710001F022730105001120181092009000800080021000800060010001000140009001300080 

1480710001F02273010600112018109200700050007001000070C070008000800080005000500030 

14807100niF022730107001120181092003000300040004000400050007000800060005000300030 

1480710001F022730108001120181092G0600U800080006000800060005000600060006000800130 

1480710001F022730109001120181092010000800060007001000090007001600110018001500140 

1480710001F022730110001120181092014000900160032001300130019000800120010001600180 

1480710001F022730111001120181092022001300280033000800080008000800220027002600170 

1480710001F022730112001120181092006000800100011000700080007000700160019001400150 

1480710001F022730113001120181092009000600070008001000100008001300260028001400160 

1480710001F022730114001120181092012001200170004000500040004000600060007000600080 

1480710001F022730115001120181092005000500120008000800100008000800100012001000080 

1480710001F022730116001120181092009000800090013002600100007000900110022000900090 

1480710001F022730117001120181092007000800080009000700070006000700070010001300090 

1480710001F022730118001120181092010001200120016001600100005000700060010000600070 

1480710001F0227301190011201810920050006000700080007000800C/000200030005000300160 

1480710001F022730120001120181092004000500070008000400070006000800080012000600050 

1480710001F022730121001120181092006000600060007000700030003000200040011001400060 

1480710001F022730122001120181092006000400020003000200070003000300020001000100060 

1480710001F0227301230011201810920G3000100020003000400030003000300130013002600350 


Because of the missing data in the first tx^ro records above, the data for 
three xjeeks, beginning with Wednesday, January 3, 1973, will be used in an 
attempt to forecast the readings for the fourth week, which is not shown 
above. The missing data value in the fourth record above will be replaced 
by 0070, x^hich is intermediate between the preceeding value 0060 and the 
subsequent value 0090. 

This data is plotted (and connected by straight line segments for 

calrity) in Figure 1. The first 32 characters in the records above are 

explained in the upper left corner of Figure 1, in accordance with the 

1 2 

standard SAROAD format ’ . The last 48 characters in the records above give 
the data which is shovm again in Table 1. 
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A linear trend line was fitted to the data with the following results: 

T = a + bX = .90 + .OOX = .90. 

The 24-hour centered moving averages were then calculated, are sho^m 
in Table 2 and are plotted in green on Figure 1, The fitting of a linear 
trend line by the least squares method to the 24-hour centered moving 
averages yielded: 

T^ ^ a^^ + = .90 + .OOX^ = .90. 

/s 

Table 3 shows D, the ratio of each two-hour average to its corre- 
sponding 24-hour centered moving average. The values of D ware averaged 
for each period of each day of the week, as shotm in Table 4. Each 

A 

average-D value was then applied to its corresponding trend value to 
produce the forecast shown in Table 5 and plotted in green on Figure 1. 

The time period average is plotted in Figure 2 for each day of the week 
and for the overall average. 

Next, the seven-day centered moving averages were calculated, are 
shown in Table 6 and are plotted in red on Figure 1. A linear trend 
fitted by the least squares method yielded: 

T 2 = a2 b^X^ = 1.00 + . 00 X 2 = 1.00. 

A A 

Table 7 contains DxW, the ratio of each two-hour average to its 
corresponding seven-day centered moving average. The values of D x W were 
averaged for each period of each day, as sho^-m in Table 8. Each average- 
DxW was then applied to its corresponding trend value to produce the 
forecast sho\\m in Table 9 and plotted in red on Figure 1. The daily 

A A 

average-D X W is plotted in Figure 3, 

No attempt is made to draw inferences from this one, short-term example. 
The only purpose here Is to illustrate the procedure which is to be pro- 
grammed. Only after a large number of such examples have been studied in 
detail should one attempt to draw conclusions about the data using this 


procedure. 









































































1973 Day; 

1 

2 

3 

HI 

5 

6 

7 

8 

9 

10 

11 

12 

00-02 

02-04 

04-06 

08-10 

10-12 

12-14 

14-16 

16-18 

18-20 

20-22 

22-24 

Wed Jan 3 







.8125 

1 

.6875 

.5708 

.4792 

.4500 

.4708 

Thu Jan 4 

.4833 

.4708 

.4542 

.4417 

.4417 

.4667 

.5083 

.5542 

.5958 

.6917 

.7750 

.7833 

Fri Jan 5 

.8000 

.8375 

,8958 

.9500 

.9958 

1.0292 

1,0250 

1.0042 

.9875 

.9375 

,8875 

.8875 

Sat Jan 6 

.8834 

.8667 

.8333 

,7916 

.7417 

.6875 

-6500 

.6250 

,6042 

.5667 

,5292 

.5084 

Sun Jan 7 

.4958 

,4917 

,4834 

.4750 

.4667 

,4583 

,4708 

.5042 

.5417 

,5667 

.5917 

.6125 

Mon Jan 8 

.6084 

.5917 

.5833 

.5875 

.6125 

.6750 

.7334 

.7500 

.7417 

.7375 

.7500 

.7708 

Tue Jan 9 

.7916 

.8417 

.9042 

.9750 

1.0542 

1.0875 

1.1084 

1.1292 

1.1750 

1.3208 

1.4375 

1.4667 

Wed Jan 10 

1.5333 

1,5500 

1,5208 

1,4917 

1.4625 

1.4834 

1.5333 

1.5833 

1.6500 

1.7042 

1.6875 

1.6458 

Thu Jan 11 

1.5792 

1.5333 

1.5750 

1.6875 

1.8000 

1.8375 

1.7667 

11.6792 

1.5833 

1.4166 

1.3208 

1.3167 

Fri Jan 12 

1.3125 

1.3042 

1,2750 

1.2167 

1.1333 

1.0750 

1.0792 

1.0834 

1.0625 

1.0375 

1.0375 

1.0583 

Sat Jan 13 

1.0708 

1.1000 

1.1667 

1.2458 

1.2833 

1.2875 

1,3042 

1.3417 

1.4084 

1.4333 

1.3958 

1.3500 

Sun Jan 14 

1.3083 

1.2625 

1.1500 

.9792 

.8584 

.7917 

.7292 

.6708 

.6208 

1 

.6167 

.6458 

.6833 

Mon Jan 15 

.7250 

.7500 

,7750 

.8125 

.8500 

.8667 

.8834 

.9125 

: .9125 

.9208 

1.0167 

1,0917 

Tue Jan 16 

1.0875 

1.0875 

1.0958 

1.1417 

1.1792 

1.1792 

1.1750 

1.1667 

1 1.1625 

1.1438 

1.0458 

.9542 

Wed Jan 17 

.9375 

.9250 

.9000 

.8333 

.8000 

.8167 

.8292 

.8584 

.8917 

.9375 

1.0042 

1.0542 

Thu Jan 18 

1.0625 

1.0583 

1.0542 

1 1.0500 

1.0208 

.9834 

,9542 

.9083 

.8625 

.8084 

.7375 

,6917 

Fri Jan 19 

.6708 

.6375 

,6042 

.5708 

.5375 

I ,5625 

I .5958 

.5875 

.5833 

.5833 

.5708 

.5542 

Sat Jan 20 

.5667 

.6083 

,6542 

.7042 

.7458 

.7125 

.6750 

.6875 

.6875 

.6792 

.6875 

,6833 

Sun Jan 21 

.6542 

.6167 

,5750 

.5542 

.5833 

.6208 

.6250 

.6167 

.5916 

.5583 

.5208 

.5167 

Mon Jan 22 

.5333 

.5375 

.5334 

.4833 

.3875 

.3333 

,3208 

.2958 

.2833 

.2833 

.2916 

.2833 

Tue Jan 23 

.2667 

.2667 

,3125 

,4083 

.5625 

.7875 
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Cable 2; 

: 24-hoi 

ur Centered Moving Averages 










1 

00-02 


2 

02-04 


3 

04-06 


06-08 


08-10 


Wed Jan 3 
Thu Jan 4 

— 

.6207 

.6372 

.6605 

.6792 

1.3584 

1.2856 

Fri Jan 5 

1.1250 

.9552 

,9600 

2.2105 

.8034 

.5830 

Sat Jan 6 

.7924 

.5769 

.8400 

1.2633 

.9438 

1.0182 

Sun Jan 7 

.6051 

.6101 

.8275 

.8421 

.8571 

1,0910 

Mon Jan 8 

.9862 

1.3520 

1.3715 

1.0213 

1.3061 

.8889 

Tue Jan 9 

1.2633 

.9505 

.6636 

.7179 

.9486 

.8276 

Wed Jan 10 

.9131 

.5806 

1.0521 

2,1452 

,8889 

.8764 

Thu Jan 11 

1.3931 

.8478 

1.7778 

1.9556 

.4444 

.4354 

Fri Jan 12 

-4571 

.6134 

.7843 

. 9*041 

.6177 

.7442 

Sat Jan 13 

.8405 

.5454 

,6000 

.6422 

.7792 

.7767 

Sun Jan 14 

.9172 

.9505 

1.4783 

.4085 

.5825 

.5052 

Mon Jan 15 

.6896 

.6667 

1.5484 

.9846 

.9412 

1.1538 

Tue Jan 16 

.8276 

.7356 

,8213 

1.1386 

2.2049 

.8480 

Wed Jan 17 

.7467 

.8649 

.8889 

1.0800 

.8750 

.8571 

Thu Jan 18 

.9412 

1.1339 

1,1383 

1.5238 

1.5674 

1.0169 

Fri Jan 19 

.7454 

.9412 

1.1586 

1.4015 

1.3023 

1,4222 

Sat Jan 20 

.7058 

.8220 

1.0700 

1.1360 

.5363 

.9825 

Sun Jan 21 

.9172 

.9729 

1.0435 

1.2631 

1.2001 

.4832 

Mon Jan 22 

1.1251 

.7442 

.3750 

.6207 

.5161 

2.1002 

Tue Jan 23 

1.1248 

.3750 

.6400 

.7348 

.7111 

.3810 

Average D 

.8868 

.7938 

.9850 

1.1336 

.9692 

.9138 


12-14 

14-16 

16-18 

18-20j 

20-22 

22-24 

,7385 

1.1636 

1.0512 

1.46081 

,8889 

.6372 

1.1804 

.9022 

,8392 

.7228 

.7742 

.8936 

.9756 

.9958 

1.4177 

.9600 

1,4648 

.9014 

1,2308 

1.2800 

1.3241 

.8823 

.9448 

.5901 

1.4868 

1.5867 

1.1076 

.8823 

.5070 

.4898 

.6818 

.8000 

.8090 

.8140 

1.0667 

1.6866 

,6315 

1.4169 

.9362 

1.3628 

1.0435 

.9545 

1.2392 

.5053 

.7273 

.5868 

.9481 

1.0937 

.4528 

.4764 

1.3895 

1.9060 

1.9686 

1.2911 

. 6486 

6461 

1.5059 

1,8313 

1.3494 

1.4174 

.6134 

.9689 

1.8461 

1.9535 

1.0030 

1.1852 

.5485 

.8944 

.9665 

1,1351 

.9291 

1.1708 

.9056 

.8767 

1.0959 

1.3032 

.9836 

.7328 

.5957 

,7714 

.9462 

1.9234 

.8606 

.9432 

.7236 

.8155 

.7850 

1.0667 

1.2946 

.8537 

,5240 

.7707 

.6956 

1.2370 

.8136 

1.0120 

.3357 

.3404 

.5143 

.8572 

.5256 

2.8870 

.8889 

1,1636 

1.1636 

1. 7668 

.8727 

.7317 

.4800 

.3243 

.6761 

1.9703 

2.6882 

1.1612 

.9352 

1.0142 

.7060 

.3530 

.3429 

2.1179 

.7908 

.8856 

1.0252 

1.2488 

1.0635 

1.1375 


Table 3: Ratio of two-hour Averages to 24-hour Centered Moving Averages 








































Reading and Time Period; 


1973 Day: 

1 

00-02 

2 

02-04 

3 

04-06 

1 4 

06-08 

5 

08-10 

6 

10-12 

7 

12-14 

8 

14-16 

9 

16-18 

10 

18-20 

11 

20-22 

12 

22-24 

Wed Jan 3 
Thu Jan 4 
Frl Jan 5 
Sat Jan 6 

1 

i : 

1 

i 





1 

' 

.7583 

.7506 

.7452 

.7571 

.7744 

.7857 

Sun Jan 7 

.7994 

.8071 

.8107 

.8161 

.8250 

.8613 

.8556 

.8786 

.8994 

.9322 

.9512 

.9536 

Mon Jan 8 

.9560 

.9589 

.9708 

.9941 

1.0190 

1.0369 

1.0410 

1.0393 

1.0405 

1.0358 

1.0292 

1.0298 

Tue Jan 9 

1.0292 

1.0256 

1.0250 

1.0321 

1.0387 

1.0434 

1.0488 

1.0506 

1.0512 

1.0500 

1.0506 

1.0542 

Wed Jan 10 

1.0560 

1.0590 

1,0726 

1.0970 

1.1161 

1.1292 

1.1422 

1.1530 

1.1660 

1.1738 

1,1744 

1.1744 

Thu Jan 11 

1.1720 

1.1690 

1.1678 

1.1690 

1.1720 

1.1768 

1.1792 

1.1768 

1.1777 

1,1810 

1.1821 

1.1845 

Frl Jan 12 

1.1887 

1.1917 

1.1952 

1.2012| 

1.2078 

1.2041 

1.2006 

1.2000 

1.2018 

1.2072 

1.2202 

1.2305 

Sat Jan 13 

1.2310 

1,2268 

1.2226 

1.2250 

1.2238 

1.2172 

1.2101 

1.2054 

1.2000 

1.1815 

1.1642 

1.1571 

Sun Jan 14 

1.1458 

1.1452 

1.1345 

1.1310 

1.1292 

1.1220 

1.1095 

1.1018 

1.0917 

1.0720 

1.0667 

1.0726 

Mon Jan 15 

1.0720 

1.0696 

1.0595 

1.0399 

1.0179 

1.0000 

.9934 

.9916 

.9887 

.9851 

.9833 

.9833 

Tue Jan 16 

,9803 

.9744 

.9637 

.9476 

.9327 

.9268 

.9244 

.9208 

.9202 

.9202 

.9166 

.9113 

Wed Jan 17 

.9083 

.9041 

.8905 

.8702 

.8559 

.8446 

,8345 

.8274 

.8172 

.8125 

.8155 

.8161 

Thu Jan 18 

,8149 

,8119 

.8083 

.8095 

.8167 

.8208 

.8196 

.8196 

.8131 

.8041 

.7976 

.7922 

Fri Jan 19 
Sat Jan 20 
Sun Jan 21 
Mon Jan 22 
Tue Jan 23 

.7875 

.6702 

.7815 

.6642 

,7738 

.6619 

.7625 

.6578 

.7506 

.6625 

1 

I 

. 7440 
.6881 

.7392 

.7315 

1 

1 

i 

j 

j 

.7232 

i 

,7131 

.6940 

.6768 


Table 6t Seven-da*.. Ce'-ltered Moving Averages 





1 

^ - - ^ 

2 


4 

5 

6 

7 

8 

9 

10 

11 

12 

1973 Day: 

00-02 

02-04 

04-06 

06-08 

08-10 

10-12 

12-14 

14-16 

16-18 

18-20 

20-22 

22-24 

Wed Jan 3 
Thu Jan 4 
Fri Jan 5 
Sat Jan 6 




1 

1 



1.0550 

1.0658 

1.0735 

.6604 

.6459 

. 3818 

Sun Jan 7 

.3753 

,3717 

.4934 

.4901 

.4848 

,5805 

.8181 

.9105 

.6671 

.5364 

.3154 

.3146 

Mon Jan 8 

.6276 

.8343 

.8240 

. 6036 

.7851 

.5786 

,4803 

.5773 

.5766 

.5793 

.7773 

1.2624 

Tue Jan 9 

,9716 

.7800 

.5854 

.6782' 

.9627 

.8626 

.6674 

1.5229 

1.0464 

1.7143 

1.4278 

1.3280 

Wed Jan 10 

.3258 

.8498 

1.4917 

2.9170 

1.1648 

1.1512 

1.6624 

.6938 

1.0292 

.8519 

1.3624 

1.5327 

Thu Jan 11 

1.8771 

1.1121 

2.3977 

2.8229; 

.6826 

.6798 

.6784 

.6798 

1.8680 

2,2862 

2,1995 

1.4352 

Fri Jan 12 

.5048 

.6713 

.8367 

.9158 

.5796 

.6644 

,5830 

.5833 

1.3313 

1.5739 

1.1474 

1.2190 

Sat Jan 13 

.7311 

.4891 

.5726 

.6531 

.8171 

.8216 

. 6611 

1.0785 

2.1667 

2.3699 

1.2025 

1.3828 

Sun Jan 14 

1,0473 

1.0478 

1.4984 

.3537 

.4428 

.3565 

,3605 

,5446 

.5496 

.6530 

.5625 

.7458 

Mon Jan 15 

,4664 

,4675 

1,1326 

,7693 

.7859 

1.0000 

.8053 

.8068 

1.0114 

1.2182 

1.0170 

.8136 

Tue Jan 16 

.9181 

,3210i 

1 

.9339 

1,3719 

2.7876 

1.0790 

,7572 

.9774 

1.1954 

2.3908 

,9819 

.9876 

Wed Jan 17 

.7707 

.8848 

.8984 

1.0342 

.8178; 

. 8288 

.7190 

.8460 

.8566 

1.2308 

1.5941 ' 

1.1028 

Thu Jan 18 

1. 2271 

1,4780 

1.4846 

1,9765 

1.9591 

1.2183 

.6100 

.8541 

.7379 

1.2436 

.7522 

.8836 

Fri Jan 19 ; 
Sat Jan 20 ! 
Sun Jan 21 
Mon Jan 22 
Tue Jan 23 

.6349 

.5968 

.7678 

.7528 

.9046 

1.0576 

1.0492 

1.2162 

.9326 

.6033 

1.0753 

1.0173 

.2706 

,2734 

.4148 

.7012 

i 

i 

.4323; 

2.3641 


Table 7: Ratio of Two-hour Averages to Seven-day Centered Moving Averages 









































B. THE COMPUTER PROGRM 


A computer program has been developed by Dr. Thomas E. Pritchard, 
co-principal investigator for this research project, for 

(1) calculation of a large number of measures of location, disper- 
sion, skexmess and kurtosis; 

(2) decomposition of an extended time series analysis model; and 

2 

(3) some x goodness of fit tests. 

This computer program is described, documented and illustrated by 
examples in the attached COMPUTER PROGRAM USER^S MANUAL. 


IV , RECOMMENDAT IONS 


It has become apparent that missing data values are quite common 
in air quality data. Such missing data values create complications 
in programming which must be studied in detail and assessed for possible 
corrections. It is essential that refinements or modifications be made 
in the computer program to compensate for missing data values. 

Extensive analysis of the air quality data must be performed, 
especially with regard to detection of unrecognized periodic patterns 
and of relationships among pollution levels at different sites. This 
must be done in order to see what generalizations may be made. An 
attempt must be made to introduce weather factors (primarily temperature 
and wind velocity) into the analysis, in order to search for a relation- 
ship among these weather factors and levels of various pollutants. 

The need for continuing refinements and modifications will become 
apparent as the actual analyses of production runs of the current computer 
program progress. And, even after the analyses of the results from the 
processed data have been largely completed, new insights will identify 
needs for further refinements and modifications of the then existing 
computer program packages. 
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Introduction 

The prograniming package described in this documentation gives con- 
siderable flexibility to the user who wishes to analyze air pollution 
data. The only assumptions made about the data file are that it is in 
SAROAD format and that the data for any particular pollution parameter, 
area, and site are stored on the file in ascending chronological order. 
As presently constituted the program handles only ’’deblocked” files 
(one logical record/physical record), but this restriction could be 
removed by some minor modifications. 

The program calculates a number of summary statistics, including 
the mean, variance, geometric mean, and coefficients of skewness and 
kurtosis. These statistics are calculated on a daily, monthly, or 
annual basis (average figures for time spans longer than one year may 
also be determined). Several types of frequency distributions are 
determined and goodness-of-f it tests of data to the normal, log-normal, 
gamma, and exponential distributions are made. 

Also the program determines multiple regression and correlation 
coefficients, if it is desired to analyse relationships between several 
variables. Using the concept of moving averages, seasonal and weekly 
indices are also determined. 

Finally, the program has limited plotting capabilities. One 
variable may be plotted against another, or a variable may be plotted 
against time. 

The documentation consists of five sections: 

(1) POTENTIAL PROGRAM CHANGES BY THE USER 

(2) USER INPUT CARDS 

( 3 ) SUBROUTINES 

(4) DETAILED DESCRIPTION OF PROGRAMS 

(5) APPENDIX 
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The first section describes a few changes that the user may easily 
make in the set of programs. The second section is the most immediately 
useful from the user’s point of view. It gives a detailed description 
of the data cards the user needs to furnish for the purpose of indicating 
to the programs exactly what type of output he wants. Five examples are 
given of various user inputs, and the type of output that results from each 
input request is sho'^im in the Appendix. 

The third section gives a short description of the major purpose of 
each subroutine. The fourth section presents a more detailed description 
of each subroutine. 

Before we conclude this Introduction the following definition is in 


order. This definition will be used throughout the rest of the documentation. 
Period T: Period T will be used to refer to the period for which the user 
requests data analysis. For example, if summary statistics are requested 
for the period 4/1/75 to 7/1/77, period T would refer to the time span 
from 4/1/75 to 7/1/77. 


Let X. , = observation j on day i, where i is measured from the beginning 
of period T. (For example would represent the 4th 
observation on the first day of period T) . 



= the daily mean for day i 


= the mean o 
n 



f observation j for all of period T (that is 

X. ,/n, where n is the total number of days in period t) . 


X . = the mean of observations j for one month, 

•J 


Potential Program Changes by the User 


There are two types of change which the user may wish to make in the 
program. Both types may be made by changing appropriate program cards in 
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the main routine. Thus these cards are in the first twenty cards of the 
deck. 

The first type of change regards the logical unit numbers used for 
input units and output units in the program. The variable IDEV repre- 
sents the number of the logical unit used to read in the air pollution 
data file. The statement: 


IDEV=8 

indicates that 8 is the number of the logical unit used. This statement 
may be changed to conform with the user’s needs. The variable IN repre- 
sents the number of the logical unit used to read in the user’s input 
cards, which describe what printed reports are desired (se® section USER 
INPUT CARDS for more detail about these cards). The statement: 


IN=5 

indicates that 5 is the number of the logical unit used. This statement 
may be changed to conform with the user’s needs. The variable lOUT 
represents the number of the logical unit used to print out all desired 
information. The statement 


I0UT=6 

indicates that 6 is the number of the logical unit used. This statement 
may be changed to conform with the user’s needs. 

The second type of change regards the maximum sizes of certain arrays 
used in the programming package. The variable MAXSE2 represents the 
maximum number of days for which data can be analyzed. The statement 


MAXS22=1464 
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indicates that period T cannot be over 1464 days (4 years) in length. If 
a longer period than this is desired, the statement must be changed. The 
variable MAXS31 represents the maximum number of variables (independent 
+ dependent) that can be used in the regression analysis. The statement 

MAXS31=7 

indicates that there cannot be more than 7 variables in the regression 
analysis. If more variables are desired, the statement must be changed. 

Also, if the value of MAXS31 or MAXS32 is changed the Dimension state- 
ments for the arrays DA0FM0, M00FYR, YEARTP, and DAXBAR need to be 
changed in the STATl routine. Likewise, the Dimension statements for 
the X, EX, and Y arrays need to be changed in the C0RREL routine in this 
case. 

User Input Cards 

The user needs only a few data dards to indicate exactly what data he 
wants to have analyzed, and what type of output he wants. The first data 
card will indicate which routine should be called (STATl or CORREL) , STATl 
giving most summary statistics j and CORREL giving correlation and regression 
information, and indices based on moving averages (see descriptions of these 
routines in the SUBROUTINES section of the documentation). Data cards after 
the first one X'jill indicate the information that should be retrieved from 
the air pollution data file and specifically the output that should be gen- 
erated. Both routines mentioned above can be called in the same run, so that 
there may be several sets of data cards of the type described belox^. 

First card - Either ’SMI’ or ’REG’ should be punched in columns 1-3, 
indicating calls to the STATl and CORREL routine, respectively. 
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I. If first card had ’SMl’ in columns 1-3: 

A, Second card - this indicates what information is to be retrieved 
from the air pollution data file. Its setup is very similar to 
the setup of cards in SAROAD format, to minimize confusion. 

Following is the information on that card. 

Col. 1 - a number that is 1 if the data are collected more 

frequently than every 24 hours; the number is 2 othert<rise 

Col. 2-3: state code 

Col. 4-7: area code 

Col. 8-10: site code 

Col. 15-16: The number of the first year for which data are to be 
retrieved from the pollution file. 

Col. 17-18; The number of the first month for which data are to be 
retrieved. 

Col. 19-20: The number of the first day for which data are to he 
retrieved. 

Col. 22-23: The number of the last year for which data are to be 
retrieved. 

Col. 24-25: The number of the last month for which data are to be 
retrieved. 

Col. 26-27: The number of the last day for which data are to be 
retrieved. 

For example, if we wanted to retrieve data for the period from 
April 15, 1974 to June 28, 1975 we would punch, in columns 15-27: 
740415 750628. 

The period described by columns 15-27 will be referred to as period 
T (see definition in Introduction). 

Col. 29-33: The code for the pollution parameter (for example, 

42401 for sulfur dioxide) . 

B. Third card - This card indici^tes what summary output is desired. 

If this card is left blank, all potential output from the STATl 
routine will be generated, wh.Lch includes the following: 

1, Listing of the information from the 80-column records from the 
air pollution data file (i.e., the "raw” data) for period T. 

2. Daily statistics, including the number of observations, daily 
mean, variance, coefficient of skewness, coefficient of kurtosis, 




mean of the natural logs of the observations, log variance, 
geometric mean, minimum, maximum, range of observations, 
mid-range, and median observation. 


3. Monthly statistics, which include the same statistics as the 
daily statistics, except that the average absolute deviation 
is calculated, and the median is not calculated. 


4. A frequency distribution for the Individual observations and 
a frequency distribution for the logs of the observations, for 
each month. For width of intervals, see writeup on FRQ 
routine. Also the chi-squared statistic, for attempted fits 
to the normal, log-normal, gamma, and exponential distributions, 
is printed. 


5. 


i) A frequency distribution for the observations adjusted by 

subtracting from each observation the daily mean of the 

observation, (i.e., X. .-X, ). Also a frequency distri- 

xj 1. t 


bution for the observations adjusted by subtracting X 
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from each observation is printed. The Chi-squared statistic 
for attempted fits of this adjusted data to the normal 
distribution is also printed. The frequency distributions 
calculated here are based on one month of data. 


ii) The means X ^ during the month (e.g., in the situation 


where there are 12 daily observations, the mean for all 
of the 2 a.ra. observations for the month will be printed 
out, the meaa for all of the 4 a.m. observations will be 
printed out, etc.) 


6. The annual summary data. Exactly the same statistics as the 
monthly statistics, except for the average absolute devia- 
tion, are printed. 

7, The same sorts of frequency distribution as in (4) above for 
all observations for an entire year. 


8, The mean for each day of the week, these means being based on 
all observations over period T. 


9. The same means as in (5ii) , except that they represent all 
observations for period T instead of one month’s observations 
(. '.e. , the means are the values X 
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10. The frequency distributions for the daily means (X^ ) and 

the logs of the daily means over all of period T. Also, the 
chi-squared statistics for attempted fits of these daily 
means to the normal, log-normal, gamma, and exponential 
distributions are printed. 


11. A plot of daily means (X ) on the Y-axis, with time on the 
X-axis . 

If the third card is not left blank, a punch of 1 in column 1 of 
the card will result in suppression of output (1) mentioned above; a 
punch of 1 in column 2 of the card will result in suppression of output 
(2) mentioned above, etc. For example a card that is blank except for a 
1 in column 6 would result in all data being printed, except for the 
annual summary data mentioned in (6) above. 

If a 2 is punched in column 3, the normal identification of area, 
site, pollutant, etc, , that appears in the upper left-hand corner of the 
monthly summary page will be suppressed except for the first month. The 
user would normally only punch a 2 in column 3, when he was suppressing 
the outputs mentioned in (2), (4) and (5) above. 

II. If first card had ’REG^ in columns 1-3; 

A. Second card - In columns 1-2 the number of variables to be useu 
(including the dependent variable) in the multiple regression 
and correlation should be punched. In column 5 a, "D" should be 
punched if the daily averages are to be treated as observations. 

If, however, the individual observations are to be used in cal- 
culating the regression coefficients and correlation coefficient 
a blank should be left in column 5. 

If, however, the first variable represents time (see third 
paragraph below) the daily averages will be treated as individual 
observations, regardless of whether a ”D” is punched in column 5. 

On the same card, the user can indicate some of the output 
options he desires, formally the output will include (1) a listing 
of the 80-column records from the air pollution data file for 
period T; (2) the multiple regression equation and multiple corre- 
lation coefficient; and (3) the plot of variable number two vs- 
variable number one (see below). If column 6 of card 2 contains 
'1', the 80-column listing is suppressed. If column 7 contains ’’1’ 
the plot is suppressed. If variable number one represents time 
then the seasonal and x^eekly indices x^ill be printed (if the total 
time period is at least 2 years) but the plot will be suppressed. 
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B. There should he an additional card for each variable in the 
analysis. Thus, if 06 was punched in columns 1-2 of the second 
card, there should be 6 cards of the type described below. These 
cards will be referred to as ’Variable" cards. If time is a 
variable in the regression and correlation, the card corresponding 
to it should be first. The card corresponding to the dependent 
variable should be the last of the cards of the type described 
below. 

The format of the card for each variable should be exactly the 
same as the format of the second card used by the STATl routine 
(described in (lA) above) . If columns 29-33 are left blank in the 
first "variable" card, then time is treated as the first inde- 
pendent variable. 

To recognize the flexibility in analysis that these "variable" 
cards allow, assume for a moment that there is one independent and 
one dependent variable. Thus there will be two "variable" cards. 

If these cards contain the same numbers, except for the site code 
in columns 8-10, then data for two different sites will be compared 
for the same pollutant and the same time period. If the pollutant 
code is the only number different on the two cards, then data for 
two different pollutants for the same site and time period will be 
correlated. If the time period is the only item different on the 
two cards, then data for the same pollutant and site will be 
compared for two different time periods. 

We will now consider several examples of these data cards. We will 

explain the meaning of the numbers and will provide sample output that 

results from, the requests by the data cards. 

Example 1: SMI 

1480710001 721101 730430 42401 

1 11 

The first card, with SMI in columns 1“3, indicates that the report repre- 
senting summary information for one variable is desired. 

Second card: ’1’ in column 1 indicates data are collected with frequency 
greater than once every 24 hours. 

'48' in columns 2-3 indicates Virginia. 

'0710' in columns 4-7 indicates Chesapeake 
'001' in columns 8-10 indicates Armory. 

'721101’ in columns 15-20 indicates that the beginning of the period for 


A9 


which we want suuanary statistics calculated (i.e. , period T) is Nov. 1, 

1972. 

'730430' in columns 22-27 indicates ti civi of the period is April 30, 

1973. 

’42401' in columns 29-33 indicates that we want summary statistics for 
sulfur dioxide for the period indicated. 

Third card: The '1’ in column 1 indicates that the listing of the raw 
data from the 80 column input recordt in the air pollution data file will 
be suppressed. That is, the raw data ^rom the air pollution file for 
the period 11/01/72 to 04/30/73, on which the summary statistics will be 
based, will not be listed. The '1' in column 8 indicates that the mean 
for each day of the week, for the period 11/01/72 to 04/30/73, will not 
be printed. The ’1’ in column 11 indicates that the daily means for the 
period 11/01/72 to 04/30/73 will not be plotted. All other summary out- 
put for the one-variable case will be printed. See Section 1 in Appendix 
for the output generated. 

Example 2;Now suppose that the first two cards are the same as in 
Example 1, and the third card contains: 

11111110110 

Since there is a '1' in every column except for co’eiiu 8 and 11, all 
output will be suppressed, except for the mean for each day of the week, 
and the plot of the daily means. See Section 2 in Appendix for the 
output generated. 

Example 3:Let us consider a final example of the use of the STATl routine. 
Assume that the third card contains: 






01211000000 


In this case the daily statistics will be suppressed, all frequency 

distributions for monthly statistics will be suppressed and item 5(ii) 

will be suppressed. See Section 3 in Appendix for the output generated. 

Example 4: Now suppose that four data cards contain the following 

information: 

REG 
02 D 

1480710001 721101 730430 42401 

1480710001 721103 730502 42401 

The first card has *REG’ in columns 1-3, which is a request for correla- 
tion/regression analysis. The '02' in columns 1-2 of the second data 
card indicates that there are to be two variables in the correlation and 

regression analysis. The 'D' in column 5 indicates that if the frequency 

of observations is greater than once per day, the daily averages will be 
calculated and used in place of the individual observations for the "raw 
data" in the correlation/regression. The information on the third card 
is exactly the same data as is on the second card in Example 1, and has 
exactly the same meaning. Note that the information on the fourth card 
is the same as is on the third card, except that the beginning and ending 
dates have been changed slightly (both have been lagged or shifted to the 
right two days). The effect is that the sulfur dioxide level for 11/01/72 
will be compared with the level for 11/03/72, the level for 11/02/72 will 
be compared with the 11/04 level,.,., the level for 4/30/73 will be com- 
pared with the level for 5/02/73, ^he third card represents variable num- 
ber 1; the fourth card, variable number 2. The output that will result is 
indicated in Section 4 of the Appendix. 

Example 5: Suppose that four data cards contain the following information 



All 


REG 
02 D 

1480710001 721101 730430 

1480710001 721101 730430 42401 

The first two cards are the same as in Example 4* Note, however, that 
the beginning and ending dates on card 4 are different than the dates 
on card 4 in Example 4. These dates are the same as they are on card 3. 
Most importantly, note that columns 29—33 on card 3, which normally contain 
the pollution parameter code, have been left blank. The blank for the 
parameter code means that the first variable will represent time. The 
data for the second variable, which is represented by card number 4, will 
automatically be aggregated into daily averages, and the daily averages 
will be treated as individual observations in the regression of this 
variable against time. 

The regression equation and multiple correlation coefficient will 
be printed, as in Example 4. Although ordinarily the seasonal and weekly 
indices will be printed when the first variable represents time, the cal- 
culation and printout of these indices are suppressed in this case. This 
is due to the fact that the total length of period T is less than 2 years. 
(Two years is the minimum period length required in order for seasonal 
and weekly indices to be calculated) . 

On both examples 4 and 5, columns 6 and 7 on the second card are 
left blank. The fact that column 6 is blank means that the data from the 
80-column input records of the air pollution data file will be listed. 

The fact that column 7 is blank in example 4 tieans that a plot of the 
second variable (on the Y-axis) vs, the first variable will be made. 
Hotjever, since the first variable in example 5 represents time, no plot 
will be generated in example 5, even though column 7 is blank. (Output 
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option 11 of the summary statistics output allows for plotting of daily 
means of a pollution parameter, with time represented on the X-axis). 

Important note: If the air pollution data is also read in via card 

file, then the user input cards described under (I) or (II) above should 
precede the air pollution data. If summary statistics calculations (as 
described under I) and regression analysis (as described under II) are to 
be done in the same run, then the cards requesting summary statistics 
should immediately precede the air pollution data to which the request 
is applicable. Likewise, the cards requesting regression should immediately 
precede the air pollution data to which the request is applicable. 


Subroutines 


1. Main program 

2. STATl 

3. DAYS 

4. FINDDA 

5. CORREL 

6. FRQ 

7. GOOFIT 

8. PL0T 

9. GMIFIT 
10* INV 

Main routine - the main routine is a very short one that Indicates what 
numbers are to represent the tt<ro input units and one 
output unit, the sizes of some of the major arrays, and 
x^hat type of analysis is to be performed on the data. If 
the data card to indicate the type of analysis has ’SMI’ 
in its 1st three columns the STATl routine will be called 
by the main program (see STATl description to indicate the 
function of the STATl routine). If the data card to 
indicate the type of analysis has 'COR* or ’REG’ in its first 
three columns the CORREL routine will be called by the main 
program (see CORREL description) . See the section ’’USER 
INPUT CARDS” to find out more about data card mentioned 
above . 


STATl - The STATl routine is used to do most of the major summary calcula- 
tions on the pollution data. It provides calculations of; 

A) daily summary data 

B) monthly summary data 

C) annual summary data 
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D) frequency distributions for: individual observations; logs of 
individual observations; adjusted observations; daily means 

E) goodness of fit to nomal, gamma, and exponential distributions 
for frequency distributions mentioned in (D) , 

F) averages for the same time each day during period T. 

G) averages for each day of the week during period T. 

It also will provide limited plotting capabilities. (See section 
"USER INPUT CARDS" and Examples 1-3 for more complete description 
of each of the items above) . 

The STATl routine calls the DAYS ,FINDDA,FRQ,GOOFIT, PLOT, and GAMFIT 
routines. For further information about these, see the paragraphs 
on each of these routines. 

DAYS - The DAYS routine receives a month number, day, and year through its 
argument list (variables IM,JD, and KY, respectively) and derives 
a value for a variable (L) which represents the days elapsed in the 
20th century. This is so that consecutive days may be recognized 
as such when filling up arrays (for example 03/01/77 is recognized as 
being 1 day later than 2/28/77. Also days of the week (Monday, 
Tuesday, etc.) may be properly ascertained when necessary (see item 
G under description of STATl routine). 

FINDDA - The FINDDA routine is referenced by both the STATl and CORREL 
routines. It is used to read the file of air pollution data, 
retrieving information from the file as requested in the STATl 
or CORREL routine. 

CORREL - The CORREL routine is used to do any correlation/regression 
calculations requested (see write-up on main routine and the 
section "USER INPUT CARDS" to see how to request correlation/ 
regression) . If the independent variable in the regression is 
Time, then moving averages are also calculated on a "seasonal" 
basis (in reality an index is determined for each day of the 
year) and on a weekly basis (i.e. , an index value is calculated 
for each day of the week). The CORREL routine calls the DAYS, 
FINDDA, INV and PL0T routines. 

FRQ - The FRQ routine will take a given set of observations, along with 
their mean and variance, and construct a frequency distribution of 
these observations. As set up now, IS intervals are used, the 
first interval representing all observations more than 2 standard 
deviations belox^ the mean. The second through 17th intervals each 
have width 1/4S (e.g. the second interval represents all observations 
between 1.75 and 2 standard deviations below the mean, etc.) The 
18th interval represents all observations greater than 2 standard 
deviations above the mean. Routine is called by the STATl routine. 


GOOFIT - The GOOFIT routine accepts the frequency distribution generated 
by the FRQ routine and calculates the chi-squared value for a 
test for goodness-of-fit of the observations to the normal dis- 
tribution* A goodness-of-fit test for the logs of the observations 
is also done. Routine is called by the STATl routine. 

PL0T - The PL0T routine will generate plots (via the line printer) of one 
variable vs. another. Routine is called by STATl and CORREL. 

GAMFIT - The GAtlFIT routine accepts the frequency distribution generated 
by the FRQ routine and calculates the chi-squared value for a 
test for goondess-of-fit of the observations to the Gamma dis- 
tribution and to the exponential distribution. Numerical inte- 
gration by Simpson’s rule is used in order to evaluate probabil- 
ities for the Gamma distribution. 

INV - The INV routine handles matrix inversion. It is called by the CORREL 
routine to invert a matrix, in the process of finding the multiple 
regression coefficients. 





DETAILED DESCRIPTION OF PROGRAMS 


Note; The meanings of a few of the more important variables are indi- 
cated by comment cards inserted in the program deck. However^ the 
general flow of the programs is indicated by the following section, 
SUBROUTINE STATl 


Cards 

710- 

720 

This statement reads the state, area, site, pollution 
parameter, and time period (period T) for which analysis 
is requested. 

Card 

730 

The It^R^PT array indicates which of the summary output will 
be printed, and which will be suppressed. 

Card 

790 

JT0TDY will be the total span, in days, of the period for 
which data are to be analyzed (period T) . 

Card 

900 

Call EINDDA to read some records from air pollution data 
file. 

Cards 

980- 

990 

If NUMDAY=0, no data for period T was found, so error 
message is printed, and execution is terminated. 

Cards 

1000- 

1040 

lOBS is the number of daily observations. 

Cards 

1620- 

3610 

After all necessary data have been read for year, or data 
have been read for the last day of period T, most of the 
statistics are calculated. Each pass through the DO loop 
250 represents one month; each pass through loop 240 
represents 1 day of the given month. 

Cards 

1320- 

1470 

Do loops 40, 44, and 48 determine area, site, and pollution 
parameter names that correspond to codes for the data that 
are being analyzed. 

Card 

1790 

N0DAYS is the number of days in given month. 

Cards 

1890- 

1960 

Each pass through Do-loop 230 corresponds to one observa- 
tion for the given day. 

Cards 

2090- 

2160 

Do loop 30 sorts the daily observations. 

Cards 

2180- 

2210 

AMED, RANGE, and AMIDRA are median, range, and raid-range 
of daily observations, respectively. 

Cards 

2220- 

2510 

Calculate daily mean, variance, log mean, log variance, 
geometric mean, skewness, and kurtosis. Also add to running 


totals for monthly statistics. 


Cards 2720 
2740 

Cards 2750' 
2900 

Cards 3100' 
3210 


Cards 3260' 
3460 


Cards 3620' 
3760 

Cards 3860' 
4050 


Card 4210 


Cards 4240- 
4490 


Cards 4670 
4940 


Cards 4960- 
5050 

SUBROUTINE 

Cards 5770- 
5820 


Card 5840 


Do-loop 60 calculates the average absolute deviation from 
the mean for the month. 

Calculate monthly summary statistics. 

^ 1 

Do loop 98 calculates the means X.j for the observations 
for the month (output item 5(ii) under Section I of 
"USER INPUT CARDS"). 

Do loop 101 adjusts individual observations for the month 
by subtracting the daily means on one pass and by subtract- 
ing the means X.^ on the second pass. 

Next 15 statements finish calculation of yearly summary 
statistics. 

Do loop 145 calculates the frequency distribution of all 
observations for the year. Also goodness-of-fit tests 
are done. 

Statements beyond this point are only executed after all 
necessary data for period T have been read. Various 
summary figures for the entire period are computed. 

Statements doim through statement 170 are to calculate and 
print out average pollution for each day of the week, for 
all of period T (these averages are output item I, B(8) in 
the description of "USER INPUT CARDS". 

The frequency distributions for the daily means and the 
logs of the daily means, for all of period T are calculated. 
Also goodness-of ~f it tests are done. The printout for this 
is mentioned in Item 1, B(10), under the description of 
"USER INPUT CARDS". 

A plot of the daily means vs. time (month, day, and year) 
is done for all of period T, 

FINDDA 

These two read statements are used to read records from the 
air pollution data file. The first read statement is used 
for records which contain a "1" in column 1 (frequency of 
observations greater than once/ilay), T’k second read 
statement is used for records that contain "2" in column 1. 

If the CORREL routine called the FINDDA routine, the rest 
of FINDDA is skipped. The checking (to see if the record 
just read corresponds to information requested; is done in 
the CORREL routine. 
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Curd* 5880 
5980 


Card 6000 


SUBROUTINE 
Card 6430 

Card 6500 
Card 6540 

Cards 6610' 
6620 


Card 6700 


Cards 6720- 
6740 

Cards 6790' 
7370 


Cards 7120 
7200 


If SAROAD data contains ”2" in column 1, there may be data 
for up to 4 pollutants on same record. Record just read is 
checked to see if pollutant for which analysi* was requested 
matches any of the (up to 4) pollutants mentioned on SAROAD 
record. 

This point is only reached if pollution parameter, site, 
and area on SAROAD record match those for which analysis 
was requested. If the date on the SAROAD record also is 
contained in period T, the RETURN statement xjill be executed, 
and data from the SAROAD record will be passed back to the 
calling routine. 

CORREL 

TEND (I) is 0 until pollution data for variable I have been 
read from the air pollution file for a date beyond the last 
day of period T for variable I. When lEND(I) is non-zero 
for all I, there is no further need to read the file. 

The X(I,J) array contains the various observations for each 
variable (I) over the time span of period T. 

NOVAR is the total number of variables. AGGREG is the 
aggregation level for the data (if AGGREG is 'D' all 
observations for one day will be averaged together and the 
daily average will be treated as one observation) , The 
array IWROPT will control the output that will be printed 
(see writeup USER INPUT CARDS). 

The Read statement reads the parameters that define each 
variable in the regression problem. Submission of a card 
that is blank (except for the indicated beginning and ending 
of the period for which data is to be analyzed) for the 
first variable will cause the first variable to be treated 
as time. 

If first variable represents time, data for dependent 
variable will be aggregated on daily basis (daily averages 
will be treated as individual observations) . 

A check is made to see if the arrays have been dimensioned 
large enough to allow for the time span (in days) of period T. 

This Do- loop checks to see if the record just read by FINDDA 
corresponds to a parameter requested, in the time period 
requested. It also fills in the X array. If the number of 
individual observations for each variable over the time period 
for which data are to be analyzed exceed the value of MAXSS2, 
all observations on one day are averaged (the daily average 
then being treated as one observation) . 

The nine statements aggregate the data for one day into one 
daily average, which is then treated as one observation. 


Cards 7220- 
7310 


Cards 7380- 
7400 


Cards 7420- 
7430 


Cards 7470- 
7570 


Cards 7720- 
7790 


Cards 7830- 
7970 


Cards 8000- 
8110 


Cards 8120- 
8130 


Cards 8140- 
8350 

Cards 8360- 
8440 


Cards 8470- 
8640 


This Do-loop is executed if the data are not to be aggregated 
into daily averages, or if the frequency of observations is 
one/day or less. 

This Do- loop checks to see whether the most recent record 
read from the tape corresponds to a date that is after the 
last date for which analysis is requested for any variable. 

This Do-loop fills up the first row in the X array 
(representing the first variable's values) when the first 
variable is time. One unit represents 30 days. Since 
X( 1,1+1) represents one day later than X(1,I) the difference 
between X( 1,1+1) and X(1,I) is 1/30. 

These Do- loops calculate the sum of the observations for 
each variable and the sum of the cross-products of the obser- 
vations. If an observation is missing for one variable, the 
corresponding observation for the other variables is discarded. 

The routine INV is called to invert the matrix of sums and 
cross-products in order to determine the coefficients of the 
regression equation and the multiple correlation coefficient. 
The coefficients of the regression equation are returned as 
elements of column NOVAR+1 of the A array. The variable R 
represents the multiple correlation coefficient. 

Do- loops 81, 84, and 88 determine area, site, and pollution 
parameter names that correspond to the codes for the data 
that are being analyzed. 

The set of Write statements in Do-loop 70 and immediately 
follox^ing that write out a verbal identification for each 
variable and then the regression equation and correlation 
coefficient. 

If the first variable does not represent time or the length 
of period T is less than two years, the seasonal and cyclical 
indices are not calculated (the "GO TO 120" is executed). 

Before doing any time series analysis, data gaps are 
"filled in" by the procedures in these statements. 

This group of statements calculates an index number for 
every day, beginning with the 183rd day of period T and 
ending with the 182nd day prior to the end of period T 
(i.e., MAXPOS-182, where MAXPOS is the total number of 
days in the period). The index is calculated by dividing 
the daily mean by the moving average centered on that day. 

These statements will check to see x^hether any day of the 
year has had more than one index calculated above. If any 
day has had two or more index numbers calculated for it 
(which x^ill be the case if the total time span for which 
analysis is requested exceeds 2 years) the index numbers 
will be averaged to derive one index for that day. 


Cards 8650- 
8840 


The statements down through statement 300 will calculate 
weekly indices in the same general manner that seasonal 
indices were calculated. 


Cards 8860- 
8900 


D 0 “loop 340 will print out the 7 weekly index numbers, 
identified by day of the week. 


Cards 8920- 
9050 


Do-loop 330 will print out the 365 "seasonal” index 
figures, which will be identified by the number of the 
day of the year. 


Cards 9070- The rest of this subroutine is for the purpose of 

9290 plotting values of variable 2 against values of variable 1, 
if such a plot is requested by the user’s input cards. 

Do loop 350 sorts the pairs of observations so that the 
values of variable 1 are in ascending order. Variable 1 
will represent the X-axis. 


SUBROUTINE FRQ 


Cards 9740- This routine calculates the frequency distribution for 

9750 observations and logs of observations. If ILOG is 0, the 
frequency distribution is based on the observations. If 
ILOG is 1, the frequency distribution is for the logs of 
the observations. 


Cards 9810- Do-3oop 80 assigns each observation (or log of observation) 
9870 to the appropriate interval (i.e. one of the class counts 
is increased by 1). There are 18 classes, or intervals, 
the width of all but the first and last interval being 1/4 
standard deviation- The first and last intervals represent 
all observations more than 2 standard deviations below, and 
more than 2 standard deviations above, the mean, respectively 


Cards 9910- 
9950 


Do-loop 96 prints out the lower and upper bound for each 
interval, and the number of observations falling in each 
interval. 


SUBROUTINE GOOF IT 

Cards 10150- Do-loop 20 computes the Chi-square statistic for goodness- 
10250 of-fit of the data to a normal distribution. For each of 

the 13 intervals or classes, the expected Frequency (EXPFR) 
of observations in that class is compared with the actual 
frequency of observations (FREQ) in that class. If either 
the actual or expected frequency is less than 5 for any class 
however, that class is combined with the next class. These 
two classes are treated as one, for the purpose of calcu- 
lating the Chi-squared statistic. The number of degrees of 
freedom is represented by the variable IDF. 


A 20 


SUBROUTINE PLOT 

Cards 10560- Do-loop 20 computes 10 values of the YM array, which will 
10570 be printed along the Y-axis to identify the scaling on the 
Y-axis. The Y-axis will represent the values of variable 
number 2. 

Cards 10610- The seven statements here regulate the number of lines the 
10670 printer carriage should be moved down before the 1st point 
(representing (X^,Y^)) is plotted. 

Cards 10680- If the 1st value of variable number 2 (that is, Y ) is 

10710 recorded (i.e,,Yj^ is not 9999.) then the 1st point xvLll be 
plotted by statement 70. Otherwise nothing will be plotted. 

Cards 10720- Do-loop 90 will take care of the plotting of all points 
10910 after X^,Y^^. 

Cards 10730- These four statements determine the number of lines to ad- 
10760 vance the printer carriage before plotting the next point. 
The total number of lines that the carriage has been 
advanced is also updated. 

Cards 10780- These five statements move the printer carriage down a line 
10820 at a time until it has been moved down a distance equal to 

the value of MOVE (which represents the number of lines down 

the X-axls between the values of X. and X,,.). 

1 1+1 

Cards 10840- One of these two write statements will be executed to print 
10850 out values of X^^ (the 1st variable) to identify the scaling 
on the X-axis. the 1st variable is time, the month, day, 

and year are printed for the scale values on the X-axis. 

Card 10900 If the ith value of variable number 2 is recorded (i.e. , is 
not 9999.) then point i will be plotted by statement 150. 

SUBROUTINE GAMFIT 

Cards 11090- The statements do™ through statement 20 are to calculate 
11610 the Chi-squared statistic for a fit to the Gamma distribu- 
tion. Eighteen classes (intervals) are used. 

Cards 11140- These 3 statements calculate the estimated parameters B and 
11160 C of the Gamma distribution, as functions of the mean and 
variance. 


Cards 11170- The statements do™ through statement 15 calculate the 
11280 constant part of the Gamma density function 
(i.e., cb+i/r(b+D) . 


Cards 11330- One pass is made through Do-loop 20 for each of the classes 
11610 or intervals. The expected frequency of observations in 
each interval is computed and compared with the actual 
number of observations in each interval, in order to cal- 
culate the Chi-squared statistic. If either the actual or 
expected frequency is less than 5 for any class, however, 
that class is combined with the next class. These two 
classes are treated as one, for the purpose of calculating 
the Chi-squared statistic. The number of degrees of 
freedom is represented by the variable IDF. 

Cards 11370- The expected number of observations in each class is 

11490 evaluated by first finding the probability of an observation 
falling in that class. In order to find this probability 
each interval or class is divided into ten subintervals. 

Then Simpson’s rule is used to numerically integrate the 
Gamma density function over the interval. 

Cards 11650- The rest of this subroutine is used to calculate the Chi- 

11870 squared statistic for a fit to the exponential distribution. 
Eighteen intervals are used and the procedure is basically 
the same as for the Gamma distribution. However, numerical 
integration is not necessary, since a closed-form integral 
exists for the exponential density function. Thus, there 
is no need to divide each interval or class into 10 
sub- intervals . 
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Al» >»OLLUIIoN analysis 
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AIR PULLUUoN AuALYbls 
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Alir PiM-LUTION ANAt-TSlS 
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